%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [images,labels] = loaddigit(imageFile,labelFile,id)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LABELS
% tic;
fid = fopen(labelFile,'r');
labelHeaderSize = 8;
header = fread(fid,labelHeaderSize,'uchar');
labels = fread(fid,inf,'uchar');
fclose(fid);

% SUBSET?
if (nargin==3)
    idx = find(labels==id);
    labels = labels(idx);
    %  fprintf(1,'Loading %d handwritten %ds\n',length(labels),id);
else
    idx = [1:length(labels)]';
    %  fprintf(1,'Loading %d handwritten digits\n',length(labels));
end;

% IMAGES
side = 28;
height = 28;
width = 28;
nPixel = side*side;
nImage = length(labels);
images = zeros(nImage,height*width);

% SLURP
colormap hot;
fid = fopen(imageFile,'r','b');
imageHeaderSize = 16;
header = fread(fid,imageHeaderSize,'uchar');
diffIdx = [idx(1)-1; diff(idx)-1];
for ii=1:length(labels)
    slurp = fread(fid,nPixel*diffIdx(ii),'uchar');
    img1 = fread(fid,nPixel,'uchar');
    img1 = reshape(img1,side,side)'/255.0;
    img2 = deskew(img1);
    %   img3 = crop(img2,height,width);
    if (labels(ii) == 6)
        subplot(1,3,1); imagesc(img1); axis equal tight; colormap hot;
        subplot(1,3,2); imagesc(img2); axis equal tight; colormap hot;
        %subplot(1,3,3); imagesc(img3); axis equal tight; colormap hot;
        drawnow;
        pause
    end
    images(ii,:) = img2(:);
end;
fclose(fid);
% toc;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESLANT

function out = deskew(in)

% GRID
[nrow ncol] = size(in);
[x,y] = meshgrid(1:ncol,nrow:-1:1);
x = x(:);
y = y(:);

% TRANSLATE TO ORIGIN
m = in(:)/sum(in(:));
ux = m'*x;
uy = m'*y;
x = x-ux;
y = y-uy;

% MOMENTS OF INERTIA
ixx = m'*(x.*x);
iyy = m'*(y.*y);
ixy = m'*(x.*y);

% DESLANT
dx = ux+(x+ixy*y/iyy);
dy = uy+y;
y = uy+y;
x = ux+x;

% ROUND
dx = round(dx);
dy = round(dy);

% REIMAGE
out = zeros(size(in));
ii = find(dx>=1 & dx<=ncol);
jj = sub2ind(size(out),dy(ii),dx(ii));
kk = sub2ind(size(in),y(ii),x(ii));
out(jj) = in(kk);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CROP

function out = crop(in,height,width)

%in = interp2(in,3);
in = in.*(in>0.3);

[h,w] = size(in);
t = h-nnz(cumsum(sum(in,2)))+1;
b = nnz(cumsum(flipud(sum(in,2))));

l = w-nnz(cumsum(sum(in,1)))+1;
r = nnz(cumsum(fliplr(sum(in,1))));

[X,Y] = meshgrid(1:width,1:height);
X = (r-l)/(width-1)*(X-1)+l;
Y = (b-t)/(height-1)*(Y-1)+t;
out = interp2(1:w,1:h,in,X,Y,'linear');
out = out/sum(sum(out));
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


